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ABSTRACT 

Precise radial velocity measurements have led to the discovery of ~ 100 extrasolar 
planetary systems. We investigate the uncertainty in the orbital solutions that have 
been fit to these observations. Understanding these uncertainties will become more and 
more important as the discovery space for extrasolar planets shifts to longer and longer 
periods. While detections of short period planets can be rapidly refined, planets with 
long orbital periods will require observations spanning decades to constrain the orbital 
parameters precisely. Already in some cases, multiple distinct orbital solutions provide 
similarly good fits, particularly in multiple planet systems. We present a method for 
quantifying the uncertainties in orbital fits and addressing specific questions directly 
from the observational data rather than relying on best fit orbital solutions. This 
Markov chain Monte Carlo (MCMC) technique has the advantage that it is well suited 
for the high dimensional parameter spaces necessary for the multiple planet systems. We 
apply the MCMC technique to several extrasolar planetary systems, assessing the un- 
certainties in orbital elements for several systems. Our MCMC simulations demonstrate 
that for some systems there are strong correlations between orbital parameters and/or 
significant non-Gaussianities in parameter distributions, even though the measurement 
errors are nearly Gaussian. Once these effects are considered the actual uncertainties 
in orbital elements can be significantly larger or smaller than the published uncertain- 
ties. We also present simple applications of our methods such as predicting the times 
of possible transits for GJ 876. 
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1. Introduction 

Recent detections of planets around other stars have spurred a wide range of research on planet 
formation and planetary system evolution. The first planet discovered around a solar type star, 51 
Pegasi b, was in a surprisingly short period orbit (Mayor & Queloz 1995). Other early planets such 
as 70 Virginis b revealed surprisingly large orbital eccentricities (Marcy & Butler 1996). Multiple 
planet systems have revealed intricate dynamical interactions such as the resonances in Upsilon 
Andromedae (Butler et al. 1999) and GJ876 (Marcy et al. 2001b). 

The future of radial velocity planet searches promises to be exciting. Ongoing large surveys 
including a broad array of nearby main-sequence stars will continue to increase the number of 
known extrasolar planets. Refinements in the radial velocity technique should continue to improve 
measurement precision, permitting the detection of planets with smaller masses. The increasing 
time span of precision observations will permit the discovery of planets with larger orbital periods. 
Several stars presently known to harbor one planet are expected to reveal additional planets in 
longer period orbits (Fischer et al. 2001). 

These advances will also bring new challenges. The 51 Pegasi-like planets could be rapidly 
confirmed by independent observers (Mayor & Queloz 1995; Marcy et al. 1997) due to their short 
orbital periods and large velocity amplitudes. However, for a given signal-to-noise ratio, the time 
required to obtain observations to confirm or refute a possible planetary candidate will typically 
scale with the orbital period of the planet. Already, one known planet (55 Cancri d) has an orbital 
period of over 14 years (Marcy et al. 2002). 

It is also more difficult to obtain precise orbital elements for planets with large orbital periods. 
While early planet candidates were routinely observed for multiple periods before publication, 
recently planet candidates have been published when observations span only a single orbital period. 
Thus, it will become increasingly easy to over-interpret early orbital determinations. A further 
challenge is that the growing precision and time span of radial- velocity measurements are expected 
to significantly increase the number of stars known to harbor multiple planets. Fitting multiple- 
planet systems requires many more free parameters, so that observational data may not tightly 
constrain all the orbital parameters or even distinguish between multiple possible orbital solutions, 
especially in the presence of significant noise and sparsely sampled data. 

These trends imply that it will become increasingly important to understand the uncertainties 
in orbital elements and other parameters derived from such observations, and thus we must use the 
best possible statistical tools to analyze radial velocity data. 

In this paper we introduce a Bayesian analysis for constraining the orbital parameters of 
extrasolar planets with radial velocity observations. The Bayesian framework considers a joint 
probability distribution function for both the observed data (d) and model parameters which can 
not be directly observed (x). This joint probability, p(d,x), can be expressed as the product of 
the probability of the observables given the model parameters, p(d\x), and a prior probability 
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distribution function, p(x), which is based on previous knowledge of the model parameters. Bayes's 
theorem allows one to compute a posterior probability density function, p(x\d), which incorporates 
the knowledge gained by the observations d. That is 

P (x\d) = = gggofe (1) 

f p(d, x)p(x) dx J p(d, x)p(x) dx 

Unfortunately, the lower integral can be extremely difficult to compute, particularly when x has a 
large number of dimensions. Even after this integral is performed, most questions require additional 
integration of p(x\d) over many of the model parameters. This paper describes the application 
of Markov chain Monte Carlo (MCMC) simulation using the Metropolis-Hastings algorithm and 
the Gibbs sampler to perform the necessary integrations. This technique allows us to accurately 
characterize the posterior probability distribution function for orbital parameters based on radial 
velocity observations. 

In this paper we first describe our model for calculating radial velocities from planetary orbital 
parameters. In §3 we summarize methods for identifying the maximum likelihood orbital solution 
for a set of observed radial velocities. Then, we discuss methods of characterizing the uncertainties 
of model parameters in §4. We focus our attention on the application of Markov chain Monte Carlo 
(MCMC) simulation to estimate the uncertainties for orbital parameters in a Bayesian framework. 
In §5 we apply this technique to analyze several published extrasolar planet data sets. Finally, we 
summarize the potential of this technique as additional long period planets and multiple planet 
systems are discovered. 

2. Radial Velocity Model 

In radial velocity surveys, the velocity of the central star is precisely monitored for periodic 
variations which could be caused by orbiting companions. Each individual observation can be 
reduced to a measurement of the star's radial velocity and an estimate of the observational uncer- 
tainty based on photon statistics. Because the observations are averaged over hundreds of sections 
of the spectrum, the observational uncertainties of most current echelle based radial velocity sur- 
veys are nearly Gaussian (Butler et al. 1998). Stellar activity can also contribute to the observed 
radial velocities. However, most of the stars investigated in this paper are chromospherically quiet 
and so that the contribution of stellar "jitter" to the observed radial velocities is believed to be 
negligible. 

2.1. Single Planet Systems 

Such radial velocity observations are able to detect the acceleration of a star due to the gravi- 
tational perturbations due to an orbiting planet. We model the motion of the planet as a Keplerian 
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orbit using the parameters: orbital period (P), velocity semi-amplitude (K), eccentricity (e), argu- 
ment of periastron (o>), and mean anomaly at the specified epoch (M D ). The perturbation to the 
radial velocity (Aw*) of a star due to a planet on a Keplerian orbit is given by 

Av*{t) = K [cos(u; + T(i)) +ecos(u;)] (2) 

A planet's true anomaly is a function of time (t) and is related to the planet's eccentric 

anomaly (E(t)) via the relation 



/T(t)\ /T+7 
arctan I - I = W arctan I — — I . (3) 

The eccentric anomaly is related to the mean anomaly (M(t)) via Kepler's equation 

E(t) - e sin (E(t)) = M(t) - M = — (t - t ) (4) 
where M Q is a constant, the orbital phase at t = which is related to the time of pericenter (t Q ). 



2.2. Multiple Planets 

For a star being perturbed by multiple planets, there is no analytic expression for the exact 
radial velocity perturbation. However, in many cases the radial velocity perturbation can be well 
modeled as the sum of multiple independent Keplerian orbits. The deviations from the simple sum of 
Keplerians model can be divided into two types: short-period interactions and secular interactions. 
The magnitude of the short-period interactions (deviations from independent Keplerian orbits on 
an orbital timescale) is often small compared to the magnitude of the Keplerian perturbation and 
the observational uncertainties. The secular interactions are typically modeled as changes in the 
Keplerian orbital parameters. While the secular interactions can have large effects on the observed 
radial velocities, the timescales are typically much longer than the time span of observations. Thus, 
we model the observed radial velocity of a star as 

V*, mo de\(t, j) = ^2 Au *.«( t ) + C 3 ( 5 ) 

i 

where Av* t i(t) is given by Eqn. 2 using the orbital parameters of the ith planet and there is an 
unknown constant velocity offset (Cj) which provides no information about planetary companions. 
In practice, high precision radial velocity surveys are typically calibrated such that different ob- 
servatories have different zero-point offsets, and hence it is important to allow observations from 
different observatories to have independent constant velocity offsets. 

Since the uncertainties from individual observational are expected to closely follow a normal 
distribution, To evaluate the goodness of fit for a given model, we calculate the usual x 2 (%) statistic, 

X 2 = ^2 ( u *,modci (tkJk) ~ v*,obs (tk,3k)) 2 l<j\ (6) 
k 
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where i& is the time of the kth observation, j/, is the observatory used for the kth observation, and 
o-fc is the observational uncertainty of the kth. radial velocity observation. 

3. Maximum Likelihood Estimators 

Given a set of model parameters (x), which are typically the masses and orbital parameters 
of the planets, we wish to compare the model predictions to the observations using the model and 
X 2 (%) statistic described in §2. The maximum likelihood estimate (MLE) of the model parameters 
xmle is obtained by finding the minimum x 2 (%mle) = m ^ n x X 2 - These are the orbital parameters 
typically reported as the "best-fit" model. Finding the set of parameters which minimizes x 2 {%) can 
be challenging, particularly for multiple planet systems. Often a combination of several methods 
is used to identify the "best-fit" orbital parameters. 

3.1. Periodograms 

Promising orbital periods can be recognized as sharp dips in a plot of the period versus the 
minimum \ 2 f°r that period with the phase and amplitude allowed to vary (a periodogram). There 
are well developed methods for determining the significance of periodicities identified in this way 
and estimating the false alarm probability (Home & Baliunas 1986). While a single sinusoid cannot 
reproduce the signal of an eccentric Keplerian orbit, a periodogram allows the rapid identification 
of any periodicities in observational data, without requiring a simultaneous fit for amplitude, eccen- 
tricity, argument of periastron, or other parameters. While the use of periodograms for an initial 
exploration of parameter space may not be optimal for identifying all possible Keplerian variations 
in the radial velocities (e.g. large eccentricities spread power across multiple frequencies), they do 
provide an efficient means for identifying potential orbital periods. Accurate estimates of the orbital 
period are valuable input parameters for subsequent algorithms, particularly the local minimization 
procedures. 

3.2. Local Minimization 

Once an initial estimate of the orbital parameters is available (x Q ), an iterative minimization 
algorithm such as Levenberg-Marquardt (LM) (Press et al. 1992) can be used to refine the model 
parameters by minimizing x 2 ■ The ML algorithm will identify only a single local minimum for 
a given initial guess of model parameters, x Q . Unfortunately, the x 2 {%) surface can have many 
local minima. In particular, for multiple planet systems, the many degrees of freedom and the 
rugged x 2 (%) surface can render LM minimization particularly sensitive to the initial guess model 
parameters and vulnerable to finding a local minimum of x 2 {%) f ar away from and much less 
probable than the global minimum. In practice, local minimization algorithm such as LM are most 
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useful for fine tuning the parameters of a local minima identified by global search algorithms. 

3.3. Simulated Annealing 

Simulated annealing generalizes iterative algorithms such as LM to reduce the risk of becoming 
trapped in a local minimum (Press et al. 1992). Random perturbations are applied at each iter- 
ation. The perturbations are initially large (high temperature) and are gradually reduced (lower 
temperature). Provided the temperature is reduced sufficiently slowly, simulated annealing can 
convert a local minimization algorithm into a global minimization algorithm. Unfortunately, the 
choice of the cooling curve is important. Additionally, the number of iterations required may be 
prohibitively large. Nevertheless, simulated annealing can be a valuable tool when performing 
global non-linear minimization. 

3.4. Genetic Algorithms 

A genetic algorithm (GA) is an optimization algorithm loosely based on biological evolution 
(Charbonneau 1995). GAs are much less likely to become trapped in local minima than local 
minimization algorithms such as LM. The disadvantage is that they require orders of magnitude 
more evaluations of the goodness of fit statistic. When fitting planetary orbits to a single set of 
radial velocity data, this can be merely an inconvenient delay for analytic models, but GAs can be 
prohibitively time consuming if it is necessary to search a large number of data sets. Once a GA 
has identified a minimum, a LM-type minimization algorithm can provide an efficient means of fine 
tuning the solution. GAs have been applied to the v And and GJ876 systems (Stepinski et al. 2000; 
Laughlin & Chambers 2001). In these cases the GAs have verified the minima found by previous 
authors. Additional minima have been identified in the case of v And (Stepinski et al. 2000). 

4. Algorithms for Estimating Uncertainties 

The previous discussion has focused on identifying the maximum likelihood or best-fit param- 
eter values. It is also important to characterize the uncertainty in the estimation of the parameter 
values. Here we discuss three methods of estimating these uncertainties: constant A% 2 boundaries, 
resampling, and Markov chain Monte Carlo (MCMC) methods. 

4.1. Constant Ax 2 Boundaries 

One method of estimating the uncertainty in orbital parameters fit to the observed data is to 
evaluate x 2 {%) & t points on a grid in the parameter space near xmle- Then confidence intervals 
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can be calculated by finding the boundary (B) along which x 2 is constant. The probability, P, that 
the parameters lie inside the boundary can be calculated by 

where Ax 2 (x) = x 2 {%) ~ X 2 (%mle), and the bottom integral extends over the entire parameter 
space. This approach is equivalent to assuming a uniform prior in f in a Bayesian framework. 
Unfortunately, evaluating these integrals often requires a prohibitive number of evaluations of 
X 2 (x), particularly when the number of parameters is large. For example, to fit the radial velocities 
caused by one, two, or three planets requires at least six, eleven, or sixteen model parameters, 
respectively. If only ten grid points are used in each dimension, then this technique would require 
10 6 , 10 11 , or 10 16 evaluations of x 2 {%)- In practice, the integrands are can be poorly behaved, 
requiring many more evaluations in each dimension. Thus, this technique is usually not practical 
for estimating the uncertainties in orbital parameters, especially for multiple planet systems. 

Brown (2004) has applied a similar method to the case of HD 72659. Brown evaluated x 2 
choosing random parameter values within a region of parameter space using Monte Carlo rather 
than choosing parameters along a grid. Still, the large range of parameter space made it impractical 
to explore the entire range of parameter space with a constant density. Thus, Brown initially 
sampled a wide range of parameter space and then manually identified several boxes of parameters 
space to sample at higher density. In total Brown evaluated x 2 a * hundreds of millions of points 
to sample the allowed parameter space for HD 72659. 



4.2. Refitting to Synthetic Data Sets 

One method of estimating the uncertainty in orbital parameters fit to the observed data is to 
apply the fitting technique repeatedly to many sets of simulated data. Since the observation errors 
are believed to be very nearly Gaussian and each radial velocity measurement has a corresponding 
uncertainty estimate, it is straightforward to construct simulated data sets by adding Gaussian 
random values to the actual data points. Each set of simulated data is meant to represent a 
possible set of measurement values. If the same fitting procedure is applied to the actual data and 
each simulated data set, then one can obtain the distribution of best-fit parameter values. 

One disadvantage of the refitting technique is that one must identify the best-fit orbital param- 
eters for each synthetic data set. In principle, one should apply a global minimization algorithm 
to each data set, but in practice the computational requirements often dictate that only a local 
minimization routine will be run for the synthetic data sets. Even when using a local minimization 
algorithm, the computational requirements can be a burden. For example, for a one planet system 
with six free parameters, x 2 must be calculated thirteen times for each iteration of LM for each 
synthetic data set. For obtaining confidence intervals roughly equivalent to 3-cr, a sample of at 
least ~ 10 5 synthetic data sets is necessary. Assuming that the local minimization routine requires 
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an average of eight iterations to converge, this amounts to ~ 10 7 evaluations of \ 2 • Still, the use of 
only a local minimization algorithm on the synthetic data sets would render this type of analysis 
vulnerable to underestimating the range of allowed parameters. If a global minimization algorithm 
were used, even over a small fraction of the possible parameter space, the number of \ 2 evaluations 
required would increase by orders of magnitude. 

Unfortunately, even then the distribution estimated via resampling may not reflect the full 
range of possible parameter values, particularly in cases where the \ 2 surface is significantly asym- 
metric around the minimum. To illustrate this possibility, we use the actual observations for the 
extrasolar planet around HD 72659. In Figure 1, we show the best-fit value of \ 2 when all pa- 
rameters except orbital period are allowed to vary as a function of the orbital period (solid line). 
We also show with dotted lines the same curve, but using synthetic data sets generated from the 
actual observations of HD 72659. For all the data sets, orbital periods shorter than ~ 1700d are 
strongly ruled out and x 2 increases very slowly for orbital periods greater than the best-fit period. 
However, the variation in the location of the best-fit period across the synthetic data sets does not 
reflect the fact that the very slow increase in \ 2 f° r longer period orbits allows a very large range 
of orbital periods. 

4.3. Markov chain Monte Carlo, the Metropolis-Hastings Algorithm, & the Gibbs 

Sampler 

Bayesian inference using Markov chain Monte Carlo (MCMC) simulations provides an alter- 
native method for estimating the uncertainty of fitted parameters. The MCMC method has been 
applied to several other astronomical data sets and problems, including spectra analysis (Kashyap h 
Drake 1998, van Dyk et al. 2001), star formation history (Fernandes et al. 2001, Panter, Heavens, & 
Jimenez 2003), object detection (Hobson &: McLachlan 2003), the Cepheid distance scale (Barnes 
et al. 2003), cluster weak lensing (Rodriguez 2003), the SZ effect (Marshall, Hobson, & Slosar 
2003), type la supernovae (Wang, Yun, h, Mukherjee 2004), and especially the cosmic microwave 
background (e.g., Knox, Nelson, & Skordis 2001, Verde et al. 2003). 

The goal of the MCMC method is to generate a chain (i.e. sequence) of states (i.e. sets of 
parameter values, x*j) which are sampled from a desired probability distribution (/(x)). Such a chain 
can be calculated by specifying an initial set of parameter values, xo, and a transition probability, 
p(x n+ i\x n ). The Monte Carlo aspect of MCMC simulation refers to randomness in the generation 
of each subsequent state. The Markov property specifies that the probability distribution for 
determination of x n +i can depend on x n , but not previous states. If the Markov chain is reversible, 
that is, if 

f {x)p{x\x') = f(x')p(x'\x), (8) 

aperiodic, and irreducible, then it can be proved that the Markov chain will eventually converge 
to the stationary distribution /(x) (Gilks, Richardson, & Spiegelhalter 1996). The requirement 
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that the chain be irreducible guarantees that it is possible for the chain to reach every state with 
non-zero probability from any initial state. 

MCMC offers a relatively efficient method of performing the integrations necessary for a 
Bayesian analysis (eqn. 1), if we can calculate a Markov chain whose equilibrium distribution 
is equal to the joint posterior probability density function for the model parameters given the 
observed data 

For application to radial velocity measurements and orbit determination, the observational 
errors are believed to be very nearly Gaussian with accurately estimated variances. Thus, if the 
data were generated by the model specified by x, then the probability of drawing the observed 
values, p(d\x), is proportional to ~ exp (— x 2 (x)/2). If we choose a uniform prior in x (p(x) ~ 1), 
then the posterior distribution, p(x\d), is also proportional to ~ exp (— x 2 (x)/2). This will allow 
us to construct Markov chains with an equilibrium distribution equal to the posterior distribution, 



The Metropolis-Hastings (MH) algorithm involves the generation of a trial state (x') according 
to a candidate transition probability function, q(x'\x n ), and randomly accepting the trial as the 
next state or rejecting the trial state in favor of the current state. The MH algorithm specifies an 
acceptance probability, a(x'\x), such that the transition probability 



is guaranteed to be reversible and irreducible, provided only that q(x'\x) allows transitions to all x 
for which f(x) is non-zero. The MH algorithm acceptance probability is 



if f(x)q(x'\x) > and a(x'\x n ) = 1 otherwise. Note that the MH algorithm does not require that 
the normalization of f(x) be known a priori. While the MH algorithm guarantees that the chain 
will converge to f(x), it does not specify when the chain will achieve convergence. 

In principle, the basic MH algorithm can be implemented as follows. 

1. Initialize the chain with some xo and n = 0. 

2. Generate a trial state, x', according to q(x'\x n ). 

3. Calculate x 2 (x') f° r trial state and x 2 (x n ) for the current state 

4. Determine the ratio f(x')/f(x n ) ~ exp (- [x 2 {x') - X 2 {x n )] /2). 

5. Draw a random number, u, from a uniform distribution between zero and one. 




(9) 



p(x'\x) = g(x'|x)a(x'|x) 



(10) 




(11) 
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6. If u < a(x'\x n ), as defined by Eqn. (11), then set x n+ \ 

7. Set n = n + 1. 

8. Go to step #2. 

The choice of q{x'\x) and deciding when a chain has converged are the primary practical 
complications. The MH algorithm can be optimized for a particular problem by the judicious choice 
of q{x'\x). Poor choices can lead to extremely inefficient sampling and hence slow convergence. 
The most efficient choice for q(x'\x) would be p(x'\d), the posterior probability distribution itself. 
However, this is not possible for this application, since the whole purpose of the Markov chain is 
to calculate the posterior distribution. In such cases, a common choice for q(x'\x) is a Gaussian 
distribution centered around x. Still, there remain important choices about the correlations and 
scale of the candidate transition probability distribution. We will address correlations later in this 
section, but first discuss the choice of scale. If the trial states are chosen with too large a dispersion 
then a large fraction of the trial states will be rejected, causing the chain to remain at each state for 
several trials and to converge very slowly. If the trial states are chosen with too small a dispersion 
(<7q), then the small step size will cause the chain will behave like a random walk, i.e., the number 
of steps required to for the chain to traverse a distance, L, in parameter space would scale roughly 
as ~ L 2 /a 2 . 

Monitoring the fraction of trial states that are accepted is one way to verify that the scale chosen 
for q(x'\x) is not too inefficient. Optimal values for the acceptance rate have been determined for 
some simple cases. For example, when x has only one dimension and the posterior distribution 
function is known to be Gaussian, then an acceptance rate of ~ 0.44 is optimal, among the class 
of Gaussian candidate transition probability distribution functions centered on the current value 
of x. For parameter spaces with many dimensions a similar analysis yields an optimal acceptance 
rate of ~ 0.25 (Gelman et al. 2003). 

Thus, it is common to choose a Gaussian q{x'\x) centered on x with some guess for the 
scale parameters. If, after running the Markov chain for some time, it becomes clear that the 
acceptance rate is significantly different than the desired acceptance rate, then the scale parameters 
are adjusted, the previous Markov chain is discarded, and a new Markov chain is begun. It may 
be necessary to repeat this several times to determine an acceptable set of scale parameters. Note 
that Markov chains which serve as the basis for altering the step size are not combined with the 
final Markov chain for inference. 

For multidimensional parameter spaces it is not always obvious how to change q(x'\x) so as to 
obtain the desired acceptance rate. Even if all the parameters are uncorrelated, there are multiple 
possible scale parameters. For this reason, we chose to use a special case of the Metropolis-Hastings 
algorithm, widely known as the Metropolis-Hasting algorithm within the Gibbs sampler. The Gibbs 
sampler generates a trial x' by altering only a subset of the parameters in x for each step. We 



= x . If u > a(x \x n ), then set 
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combine the Gibbs sampler with a Gaussian candidate transition function, i.e., for the parameters 
in (x) to be updated, the candidate transition probability function is 



qWJxfi) = i ex P 



(12) 



for valid x'^ (i.e., if the model dictates that x'^ be positive definite, then trial states with negative 
x'^ are rejected). We use the index \x to distinguish elements of the vector of parameters, and the 
index n to indicate the n-th step of the Markov chain. Each (3^ is a parameter which controls the 
size of the steps for the parameter indicated by ji. Thus, our acceptance probability reduces to 

^^^{.p^),} (13) 

for valid x', and a{x'\x n ) = for invalid parameter values. 

While there are several variations, we choose the parameters to be altered in the next trial 
state according to randomly generated permutations of the model parameters. If the parameters 
are chosen one at at time, then it is easy to monitor the acceptance rates for steps involving each 
parameter separately. If any of these acceptance rates differ from the desired acceptance rate, then 
the guilty scale parameters, (3^, can be identified and adjusted for calculating the next Markov 
chain more efficiently. 

The other major difficulty in applying MCMC is determining how long the Markov chain should 
be before using it for inference. In practice, one performs multiple tests which can positively identify 
chains which have not converged. The failure of such tests to demonstrate non-convergence suggests, 
but does not prove, that the chain has converged (Chen et al. 2002). In our experience, we have 
found it advisable to construct multiple Markov chains for comparison. If all the chains have 
converged, then the distributions of all quantities of interest should be similar to within statistical 
uncertainties (Gelman & Rubin 1992). 



4.4. Comparison 

Both refitting to simulated data sets and MCMC provide reasonable methods of assessing the 
uncertainty in model parameters. Refitting to simulated data estimates the distribution of the MLE 
for the simulated data, which is not necessarily the most desirable estimator. Transcribed into the 
language of Bayesian statistics, the MLE is not necessary the mean or median of the parameter's 
marginal distribution. MCMC samples from the joint posterior distribution given a prior and the 
observational data. 

While changing the prior would require completely redoing an MLE analysis, importance sam- 
pling can be applied to a realization of a Markov chain to compute confidence intervals for several 
different priors simultaneously (Press et al. 1992; Gilks et al. 1996). This could become valuable 
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when considering alternative formation scenarios (e.g. alternative distributions for initial eccentric- 
ities and inclinations) or when additional constraints become available (e.g. new observations or 
long term orbital stability studies). Note that for importance sampling to be accurate, the posterior 
distributions under the different priors must not be too dissimilar. 

Further, calculating the best fit model for each set of simulated data requires repeatedly solving 
a complex and time consuming non-linear minimization problem. Similarly, most techniques which 
are less sensitive to local minima (e.g. GAs) are less efficient than MCMC for multiple planet 
systems which have a large number of free parameters. Calculating a longer Markov chain is a 
relatively simple task and computationally more efficient. Thus, we expect the MCMC technique 
to generalize to multiple planet systems much better than refitting. 

While both methods may be practical when fitting to a small number of model parameters, 
MCMC is particularly useful for high dimensional parameter spaces. This advantage is particularly 
important for long period planets and multiple planet systems where multiple free parameters which 
can be traded off against each other to obtain similarly good fits. The radial velocity observations 
of such systems can result in large valleys in the x 2 {%) surface which permit a broad range of 
parameter values. Within these valleys, there is often a rugged x 2 {%) surface, meaning that there 
are many nearby local minima for x 2 (x). While the rugged x 2 {%) surfaces present difficulties for 
local minimization routines, the Markov chain Monte Carlo method is able to jump between these 
local minima and accurately calculate the posterior probability distribution for model parameters. 
It is important to realize that if there were multiple local minima separated by a barrier region 
of parameter space which resulted in significantly larger values of x 2 > then Markov chain Monte 
Carlo would have great difficulty sampling from the posterior distribution unless special care was 
taken to identify the well separated local minima and a candidate transition probability function 
was carefully chosen to allow transitions between the local minima. Thus, it is still important to 
conduct a single global search over the parameter space to recognize if there are distinct orbital 
solutions separated by a large x 2 (%) barrier. 

In summary, MCMC has the following advantages compared to the other methods we have 
discussed. 

• MCMC naturally allows for correlated and non-Gaussian uncertainties in fit parameters. 

• The results can be simply interpreted as the joint posterior distribution for a given prior and 
set of observational data. 

• In some cases the resulting distribution can be efficiently updated to account for additional 
observations or alternative priors, provided that they do not greatly alter the posterior dis- 
tribution. 

• Calculating the next step in a Markov chain is much faster than performing an additional 
minimization with resampled data. 
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• Computationally MCMC is more efficient than other techniques, particularly for high dimen- 
sional parameter spaces. 

5. Example Applications 

To illustrate the application of the MCMC method to fitting radial velocity data, we will 
present several example applications to the presently known extrasolar planets. For calculating the 
Markov chain we use a set of parameters, 

x = (log-P, log A, ecosw, esincj, M a ,C), (14) 

where C is the set of mean velocity offsets (one offset for each observatory). It is common practice 
to choose a "non-informative" prior that is uniform in the logarithm of a positive definite magni- 
tude (e.g., orbital period and velocity semi-amplitude), as suggested by several scaling arguments 
(Gelman et al. 2003). Additionally, we found that the use of logP and log A instead of P and A 
increased the rate of convergence for systems where the orbital period was not tightly constrained. 
Similarly, we found that the use of esina; and ecosw instead of e and uj significantly increased the 
rate of convergence for systems with small eccentricities where u> is not tightly constrained. 

The distribution of published orbital parameters is roughly consistent with being uniform in 
log.P, log A, e, u, and M Q (aside from cutoffs at small orbital periods and large planetary masses 
and not including the very short period planets for which tidal forces have presumablely circularized 
the orbits). Therefore, for the histograms and contour plots presented below, we have sampled from 
the Markov chains using importance sampling so that they correspond to priors which are flat on 
(logP, log A, e, w, M ). 

For each planetary system, we initialize multiple chains with parameter values near the pub- 
lished values of the orbital elements or the best-fit models identified with a GA. We also start chains 
with initial parameters drawn from a Gaussian distribution centered on the published values with 
standard deviations three times the published uncertainties. In cases where our Markov chains 
indicated uncertainties greater than the published values, we constructed additional chains with 
initial parameters chosen over an even wider range. 

For these simulations, we have used a Gaussian candidate transition probability function 
(eqn. 12) and a Gibbs sampler which randomly chooses to update one or two parameters at 
each step with equal probability. (When two parameters have significant correlations, perturb- 
ing two parameters in a single step can allow the Markov chain to explore parameter space more 
quickly.) All parameters other than the one or two chosen to be updated are left unchanged. The 
order of the parameters to be altered is determined by a randomly generated permutation. We 
monitor the acceptance rates for steps involving each parameter separately. If a chain had an ac- 
ceptance rate below 0.25 or above 0.55, then we adjusted the relevant /3 M and start computing a 
new chain. To reduce the dependence on the initial parameter values, we discard the initial 10000 
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states or the first 10% of the chain, whichever is longer. Since consecutive states can are often 
highly correlated, we lose little information by consider only every m-th set of parameter values 
with m greater than or equal to ten times the number of free parameters. 

We test for non-convergence of each chain by comparing the marginalized distribution of each 
parameter in the first half and full chain. We also compare the results of multiple chains run with 
multiple values of (3. This helps us recognize chains where an inappropriate step size may have 
resulted in a chain becoming trapped in a small region of parameter space for the duration of 
the chain. Additionally, we run multiple chains with widely dispersed initial conditions to verify 
that the Gelman-Rubin test statistic, is consistent with convergence. Since we initialize the 
chains with wide range of initial parameter values, it is easy to recognize that the distributions 
of parameters based on only the early portion of the chains are dependent on the values chosen 
for the initial state of the Markov chain. The Gelman-Rubin test statistics can determine if the 
Markov chains have yet to converge by comparing the variance of each parameter within each chain 
to the variance of the parameter across multiple chains. For Markov chains which converge, the 
Gelman-Rubin test statistic approaches 1.0 from above. For our results we have checked that it is 
always less than 1.2 and usually less than 1.1 (Gilks et al. 1996). The plots presented are based on 
only the Markov chain with the most desirable acceptance rates. 

These very conservative choices are computationally inefficient, but the inefficiency can be 
tolerated for the purposes of this paper. For each system, we have computed tens of Markov 
chains, including a wide range of |/?|. Each of our Markov chains typically contains 10 6 — 10 10 steps, 
depending on the number of fit parameters and the choice of \(3\. (When \(3\ results in too high or low 
an acceptance rate, convergence requires many more steps.) Ideally, one would attempt to improve 
the computational efficiency while maintaining confidence that the chains have converged. Based 
on our experience from the systems studied in this paper, we believe that one could construct only 
five Markov chains with a single value of \(3\ (chosen to give a desirable acceptance rate) for each 
system and still be reasonablely confident that the Markov chains had converged. We believe that 
further refinements to the candidate transition probability function could make MCMC even more 
efficient. Using additional changes of variables could reduce covariances between model parameters 
could speed convergence and reduce the autocorrelation of the states in each chain. We hope to 
present such optimizations in a future paper. 



5.1. Eccentricities 

Despite the nearly Gaussian uncertainties for the radial velocity data, the uncertainties for 
orbital parameters can be significantly non-Gaussian. In particular, since the orbital eccentricity 
is confined between zero and one, we might expect significant deviations from normality when 
the eccentricity is small or nearly unity. In Fig. 2 we show the distribution of eccentricities and 
arguments of periastron for HD 76700, HD 68988, and HD 4203 based on the published observations 
and uncertainties (Tinney et al. 2002; Vogt et al. 2002). There are several interesting features in 
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these plots: 

1. For HD 76700, the distributions of orbital elements calculated by resampling and MCMC are 
similar. We find the distribution of eccentricities barely includes e = (assuming a uniform 
prior in x as defined in Eqn. (14)) and excludes e < 0.017 at the 90% confidence level, in 
contrast to the published solution, e = 0±0.04. While the published orbital solution does not 
constrain u, our simulations constrain the argument of pericenter to —100° < to < 60° (90% 
confidence interval). Additionally, the distribution of u has a noticeable asymmetry. For 
planets with small eccentricities such as HD 76700, the variables esinu; and ecosu-> generally 
have smaller uncertainties and are more Gaussian than the uncertainties in e and to. 

2. For HD 68988, the eccentricity distribution is clearly separated from zero, and the hypothesis 
that e < 0.1 is strongly rejected within the single planet Keplerian model. Again, resampling 
and MCMC methods result in similar uncertainty estimates. 

3. HD 4203 has a significant eccentricity, however the published observations still permit strongly 
non-Gaussian uncertainties in e. In this case, MCMC results in an uncertainty distribution 
with a more significant high eccentricity tail and a narrower distribution in uj than is predicted 
by resampling. 

HD 80606 has a very large eccentricity (e ~ 0.93) and a very small radial separation at 
periastron (r p ~ 0.035AU). While the orbital period is very well constrained, e, and thus radial 
separation at periastron (r p ) have significant uncertainties (Fig. 3), especially in the context of tidal 
dissipation, since the rate of tidal circularization is a steep function of r p (Rasio et al. 1996). Our 
MCMC simulations show that the current observations still permit orbital solutions with pericenter 
distances as small as a stellar radius. It is interesting to note that MCMC and resampling methods 
result in significantly different distributions for r p . The high correlation between variables does slow 
the convergence of the Markov chains, but we have run several chains for an extended duration and 
none result in a significant probability of r p > 0.045AU. If additional observations were to constrain 
the pericenter distance, this system might eventually provide an interesting test for theories of tidal 
circularization and orbital decay. 

5.2. Long Period Systems 

Constraining the orbital parameters of long period systems is a significant challenge, since 
good solutions typically require observations spanning at least one orbital period. For eccentric 
systems, it is much more difficult to identify optimal times to observe the system while still during 
the first or second orbit. In Fig. 4 we show the joint distribution of the period and eccentricity 
for the system HD 72659. Note that the MCMC simulations reveal a much larger uncertainty in 
orbital period and eccentricity than suggested by resampling or by the published error bars. After 
this plot was prepared, Marcy (2003) informed us that more recent observations of HD 72659 yield 
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a best-fit of P = 5828 d and e = 0.47 (large solid point), demonstrating the potential for MCMC 
to provide more accurate uncertainty estimates than resampling. 

In Fig. 5 we show the marginalized probability distributions for P, K, e, and u for several 
of the known extrasolar planets with the longest orbital periods. While the parameters are well 
constrained for some systems, significant uncertainties in the periods and eccentricities appear 
to be common. Both the mean value and uncertainty of orbital elements determined by MCMC 
simulations can be significantly different from the values and uncertainties suggested by resampling, 
particularly when at least one parameter has a significantly skewed distribution. 

5.3. Multiple Orbital Solutions 

For some systems the observational data do not determine a unimodal orbital solution; in 
particular there can be two similarly good fits that are well separated from each other. This is 
particularly common in multiple planet systems. In Fig. 6 we show two possible orbital solutions 
for 47 UMa. One corresponds to the published orbital solution and the other corresponds to a 
solution in which the outer planet has a much longer orbital period and large eccentricity. 

We have begun conducting MCMC simulations to estimate the parameter distributions for 47 
UMa and other multiple planet systems. For multiple planet systems several orbital parameters 
can be highly correlated greatly slowing convergence. We are exploring changes of variables which 
could accelerate convergence and permit more detailed investigations of the uncertainties in the 
elements of such multiple planet systems. 

5.4. Interacting Multiple Planet Systems 

The most easily observable deviation from Keplerian motion is generally expected to be pre- 
cession of the longitude of periastron. This suggests we adopt a precessing Keplerian model 

Au*,i(t) = Ki [cos (Ti(t) + L0 O:i + Uit) + ej cos (u 0ji + toit)] , (15) 

for each planet, where u>i is the precession rate of the longitude of periapse of the ith planet and t is 
the time of the observation. Such a model incorporates the dominant observable perturbation to the 
basic Keplerian model while introducing only one additional parameter per planet. This model is 
particularly good when the perturbations excite a single dominant secular eigenmode, as in GJ876, 
in which case the values of Cj{ should be nearly equal and time- independent. Further, in this model, 
Co can be fit to radial velocity measurement without time consuming N-body integrations. Although 
radial velocity measurements constrain only msini, the precession rates can in principle determine 
the planets' inclinations and masses. Of course, this approach assumes that there are no additional 
perturbations, for example due to general relativistic effects, the quadrupole moment of the star, 
or additional companions. 
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The GJ 876 system is an example of a interacting system in which mutual perturbations 
are detectable. Modeling GJ 876 with two precessing Keplerian orbits and applying MCMC, we 
investigate the precession rates of the two orbits. In the top panel of Fig. 7, we show the probability 
distributions for to\ (dotted line), C02 (dashed line), and the instantaneous average of the two (solid 
line). In the lower panel we plot the probability distribution for the difference of the precession 
rates. In this case, the two precession rates are comparable, allowing for a secular resonance in 
which Aoj remains small (Lee & Peale 2002). Based on the published observations, we can not yet 
detect a significant difference between the precession rates of the two planets. Rates differing by 
more than ~ 10° /yr can be excluded according to this model. 



5.5. Possible Transit Times 

Radial velocity measurements can be used to identify when a transit could occur, if the system 
were to have a favorable inclination. This can greatly aid observers who can concentrate their efforts 
at the appropriate times. Correlations between best-fit orbital parameters can have a significant 
impact on the search for transits. Thus, it is important that observers understand the uncertainty 
in the predicted transit times. Particularly, for strongly interacting multiple planet systems, there 
can be significant uncertainties in the possible transit times. 

We have fit a model based on Eqn. (15) to the published radial velocity data for GJ876. Us- 
ing a genetic algorithm we found a model with precession rates of — 44.0°/yr and — 46.8°/yr for 
the inner and outer planets, respectively. These are similar to the precession rates found for the 
best-fit orbital solution obtained via self consistent n-body integrations and local minimization 
(Laughlin & Chambers 2001). Clearly, this comparison does not include important short period 
perturbations to the orbits which should be included when forecasting transits and planning ob- 
servations. Additionally, a transit search would want to target the time of ingress or egress which 
introduces additional dependences on the stellar radii and orbital parameters. Nevertheless, we 
use this model to demonstrate the uncertainty in potential transit times in an interacting system 
and the applicability of MCMC to predicting potential transit times. The actual system is most 
likely precessing at a similar rate and the unknown precession rate contributes to the uncertainty 
in predicting potential transit times. In Fig. 8 we show the probability distributions for the center 
of the transit of GJ 876 b approximately two months after the last published observation, assuming 
sini = 1. The model which includes the precession rates (solid line) allows for a significantly wider 
range of transit times than the model which assumes fixed periapses. This significantly increases 
the duration over which GJ 876 must be monitored to ensure that the ingress or egress falls within 
the window of observations. 
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6. Conclusions 

Markov chain Monte Carlo methods provide a valuable tool for matching planetary orbits to 
radial velocity data. These techniques can be particularly useful for characterizing the uncertainties 
and correlations in estimated orbital parameters for extrasolar planets. Presently, several systems 
have large uncertainties in orbital parameters, particularly the long period planets and multiple 
planet systems. Our Markov chain Monte Carlo simulations reveal that the published orbital 
solutions can significantly over- or underestimate the uncertainty in orbital parameters. In addition 
to quantifying the uncertainties of various orbital parameters, Markov chain Monte Carlo allows 
specific questions to be investigated directly from the radial velocity data itself bypassing fits to 
orbital parameters. As more multiple planet systems are discovered, such methods could become 
increasingly valuable. 
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Fig. 1. — We show the minimum value of x 2 varrying all parameters except the orbital period 
versus the orbital period for HD 72659. The solid line is for the actual observations and the dotted 
lines are based on resampled data. Note that the very slow rise in x 2 {P) f° r larger orbital periods 
allows a much greater range in the orbital period than suggested by the variation in the minimum 
of the x 2 (P) curves for resampled data. 
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Fig. 2. — We show the distribution of eccentricities (far left), longitudes of periastron (center left), 
and two related quantities, esincj (center right) and e cos to (far right), for the companions to HD 
76700 (top), HD 68988 (middle), HD 4203 (bottom). The solid lines represent the results of our 
MCMC simulations, while the dotted lines show the results of our fits to resampled data based on 
the published observations and uncertainties (Tinney et al. 2002; Vogt et al. 2002). 
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Fig. 3. — The probability distribution for the periastron distance of the companion to HD 80606. 
Here we assume a stellar mass of 1.1M & and inclinations corresponding to sini = 1. The solid line 
shows the results of MCMC simulation and the dotted line shows the results of fitting to resampled 
data. Considering the strong dependence of the tidal circularization rate on the periastron dis- 
tance, the uncertainty in the pericenter distance makes it impossible to constrain theories of tidal 
circularization based on the present observations. 



-23- 



0.8 



0.6 - 



0.4 - 



0.2 - 






1000 



10 4 



P (d) 



Fig. 4. — Here we show the probability distribution marginalized over all variables except the 
period and eccentricity of HD 72659. The solid contours show the 1, 2, and 3-sigma confidence 
intervals, defined to contain 68.3%, 95.4%, and 99.73% of the probability distribution, based on 
MCMC simulations. The dotted contours also enclose the 1, 2, and 3-sigma confidence intervals, 
but based on resampling methods. The point with error bars indicates the published orbital solution 
and uncertainties. Both sets of contours assume a uniform prior in logP and e. The large point 
without error bars indicates the best-fit orbital solution based on more recent observations made 
after the rest of this plot had been prepared (Marcy 2003). Unfortunately, there are often large 
uncertainties in the orbital elements derived for long period planets. 
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Fig. 5. — Here we present the probability distribution marginalized over all but one variable for 
the eight extrasolar planets with the longest orbital periods (we restrict ourselves to systems with 
only one known planet). The solid line shows the results of our MCMC simulations. The dotted 
line shows the results of our fits to resampled data. The area under each curve is normalized to 
unity, however the calculated distributions sometimes have tails which extend beyond the domain 
of these figures. 
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Fig. 6. — Here we present the published radial velocity data for 47 UMa along with two fits to 
the data each employing two non- interacting planets on Keplerian orbits. The solid curve is the 
published solution, while the dashed solution is for a much larger period (three times the published 
best-fit period) and eccentricity (e ~ 0.8). The \ 2 °f the alternative fit is actually less than that of 
the published fit, but the difference is not significant. While the dashed solution is not favored by 
orbital stability requirements, 47 UMa illustrates that the observations on multiple-planet systems 
can sometimes be well modeled by very different orbital solutions. 
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Fig. 7. — In the top panel we show the probability distributions for the precession rates of the 
periapses of GJ 876 b (dotted) and c (dashed). The solid line is the instantaneous average of the 
two precession rates. In the lower panel, we show the distribution for the difference in the precession 
rates. 
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Fig. 8. — Here we show the times of the center of potential transits for GJ 876 b approximately 
two months after the last published radial velocity, assuming s'mi = 1. The dotted line shows 
the prediction using a superposition of two non-interacting Keplerian orbits to model the data, 
while the solid line allows each orbit to precess at an arbitrary rate. While our model does not 
include the short period perturbations important for actually forecasting potential transits, it does 
demonstrate the uncertainty of transit times in an interacting system. 



